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Abstract 

The microcanonical statistical mechanics of a set of self-gravitating particles 
is analyzed in a mean-field approach. In order to deal with an upper bounded 
entropy functional, a softened gravitational potential is used. The softening 
is achieved by truncating to N terms an expansion of the Newtonian potential 
in spherical Bessel functions. The order is related to the softening at short 
distances. This regularization has the remarkable property that it allows for 
an exact solution of the mean field equation. It is found that for N not 
too large the absolute maximum of the entropy coincides to high accuracy 
with the solution of the Lane-Emden equation, which determines the mean 
field mass distribution for the Newtonian potential for energies larger than 
Ec ~ — 0.335GM^/i?. Below this energy a collapsing phase transition, with 
negative specific heat, takes place. The dependence of this result on the 
regularizing parameter is discussed. 

I. INTRODUCTION 

The statistical mechanics of self-gravitating systems is amazing. It has been studied 
since long ago by Antonov Lynden-Bell and Wood Thirring [Q, and Kiessling [Q, 
among others 0. One reason for the interesting and peculiar behavior of these systems is 
that they are thermodynamically unstable. The usual thermodynamical limit exists only 
for those systems which are thermodynamically stable ^ . For a system of Np classical par- 
ticles interacting via a two body potential (f){r), a sufficient condition for thermodynamical 
stability states that there must exist a positive constant Eq such that, for each configuration 
{ri, . . . , Tat }, the following inequality is obeyed 0: 
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In contrast, self-gravitating systems do not possess a proper thermodynamical limit. More- 
over, due to the short distance singularity of the gravitational potential, the entropy is not 
even well defined: it diverges for any value of the energy . To define the thermodynamics 
of these systems the potential must be regularized at short distances. This can be done in 
many different ways. Particles endowed with a hard core is one possibility In this 

case, the potential is repulsive and singular at short distances. Other popular choices are 
the so called softened potentials, which are smooth at the origin. As shown by Thirring [§, 
the thermodynamical instability is caused neither by the singularity nor by the long range 
nature of the potential, but is due to the fact that the potential is always attractive^. The 
essential common feature of these purely attractive potentials is the appearance of a phase 
transition separating a high energy homogeneous phase (HP) from a low energy collapsing 
phase (CP) [0-|TT|. The phase transition takes place in an energy interval with negative 
microcanonical specific heat. From the dynamical point of view, both phases are also dif- 
ferent: the single particle motion is superdiffusive in the CP and ballistic in the HP ]T^,IH 



The dynamics and statistics of simple low dimensional models with long range attractive 
forces has been studied in [Il2|-[T5|]. Their conclusions support the idea of a collapsing phase 



transition as in the Thirring model. If angular momentum is conserved, the situation could 



be notably altered [|l6 



As mentioned, the usual thermodynamical limit does not exist for unstable systems. To 
have well defined thermodynamics when the number of particles Np is huge, the following 
scaling must be considered: the potential energy is rescaled by 1/Np, and then the energy 
and entropy scale with Np. It has been proved for the canonical ensemble that this scaling 
reproduces mean field theory exactly in the limit Np ^ oo . This means that correlations 



among two or more particles vanish, and therefore the equilibrium state is characterized 
by a one particle density only, which minimizes the free energy functional. Although we 
are not aware of any rigorous proof, we shall assume here that the same holds for the 
microcanonical ensemble, changing minimization of the free energy by maximization of the 
entropy functional. 

If the troubles caused by the short distance singularity are ignored, it is possible to write 
down a mean field entropy functional for self-gravitating systems, which depends only on 
the particle density. This functional is not upper bounded, and, therefore, has no absolute 
maximum, a reflection of the fact that the entropy is not defined in the finite system. 
For energies larger than ~ — 0.335GM^/i? there is however a local maximum. Below 
this energy, no local maximum of the entropy exists [Q]. This fact was explained in terms 
of a transition from the homogeneous isothermal sphere behavior to the CP at E^. The 
transition produces negative specific heat, and was called the gravo-thermal catastrophe [^. 
Very recently, it has been pointed out that the low energy phase might be described by 
a spherically non-symmetric deformation of the singular solution of the isothermal Lane- 



In the case of hard core particles the potential is repulsive at short distances and the thermody- 
namical instability is actually due to the long range forces. 
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Emden equation 



The gravitational potential must be modified at short distances to make equilibrium 
statistical mechanics applicable to self-gravitating systems. As shown by Kiessling for the 
canonical ensemble @], in the limit where the classical gravitational potential is recovered the 
equilibrium state approaches a particle distribution with all particles collapsing at a single 
point. The behavior of the system will depend on the scales at which the regularization is 
effective. There might be regularized potentials which, in the mean field limit, produce the 
global maximum of their associated entropy close to the solution of the isothermal Lane- 
Emden equation for energies E > E^.. If this is the case, a collapsing transition should occur 
at some energy close to E^.. The CP is expected to be very sensitive to the details of the 
regularization at short distances and the HP almost insensitive to it. On the other hand, 
if the regularization is effective only at very short distances, the solutions of the isothermal 
Lane-Emden equation will be global maxima of the entropy only at very high energies, 
and therefore the collapsing transition will take place at some energy much larger than E,,. 
The critical energy will be higher the smaller the scale at which the regularized potential 
differs significantly from the unregularized one. A completely similar picture was rigorously 
established by Kiessling for the canonical ensemble Q]. 

In this paper we introduce a convenient new softening procedure for the regularization 
of the gravitational potential, and we investigate the consequences in the microcanonical 
thermodynamics of self-gravitating systems. The rest of the article is organized as follows: 
in Sec. II we introduce the family of potentials to be studied; in Sec. Ill we derive the mean 
field equation and its general solution in terms of a set of algebraic equations; Sec. IV is 
devoted to the discussion of the results and Sec. V to summarize the conclusions. 



II. SOFTENED POTENTIAL 



As mentioned in the introduction, the short distance singularity of the gravitational 
potential causes many troubles. What is more, for a real system such singularity is not 
physical, since at short distances new physics must be taken into account. Thus, the potential 
should be modified at short distances to avoid the singularity. In simulations of cosmological 



problems a widely used potential is the so called Plummer softened |1T9|J20|| : 



(2) 



For r ^ cr, (Q) coincides with the gravitational potential. Other softened potentials are 
those known as spline softened |^ . The equilibrium thermodynamics of systems with these 
softened potentials has been studied in [^|, and the dynamical effects of softening were 
considered in p3|. The form of the potential at short distances is arbitrary to a large extent. 



since we do not know how the new interactions modify it. 

The problems with the singularity of the gravitational potential in statistical mechan- 
ics disappear if the equilibrium distribution is modified appropriately. For instance, if the 
approach to equilibrium is coUisionless, via violent relaxation, the equilibrium state is de- 
scribed by the Lynden-Bell statistics whose one-particle distribution is of Fermi-Dirac 
type, and produces an effective repulsion at short distances 
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The regularized potential we propose, which, as we shall discuss, has several remarkable 
features making it very convenient for thermodynamical purposes, is based on the following 
identity^: 

The singularity at the origin is removed by truncating the series to a given order A^. Let us 
consider a system of Np particles confined within a sphere of radius in 3 dimensions. The 
maximum distance between two particles is 2R. Hence, our potential must represent 1/r for 
distances ^ r < 2R. Therefore, we choose the following interaction energy between two 
particles of mass m located at r and r': 

0(|r - r'l) = -—- 2 Mr - r'\/R) , (4) 
^ k=i 

where = sm.{(jjkx) / {ukx) are spherical Bessel functions of order zero and Uk = {2k — 

1)71/2. Fig. |l] displays the singular potential and the regularized potential with = 10 and 
= 20. A similar expansion has been used to introduce simple models in low dimensions 
which allow to perform numerical simulations of systems with long range attractive forces 
with CPU time growing only as the number of particles [P^^-^ . 



What is remarkable of (H) is that each term obeys the following differential equation: 



+ ^,)Mr/R) = 0, (5) 



so that the potential (H) verifies 



where 



VnM) = 0, (6) 



N 2 

^^A^ = n(v^ + (7) 

k=l ^ 

This relation will prove very useful in the mean field analysis of the next section. 

III. MEAN FIELD ANALYSIS 

It is well known that long range forces suppress fluctuations, and thus in these cases a 
mean field analysis is accurate or even exact. We expect that, dealing with an unstable 
system in the scaling regime described in the introduction, the description of the thermody- 
namical state in terms of a one particle density, neglecting two or more particle correlations, 



gives the essential physical behavior [17|. We will derive in this section the form of the 
mean field equation and of the corresponding thermodynamic quantities for a system whose 
dynamics is governed by a potential of the form (^. 



^Equation (^) follows immediately from the sine series expansion of the constant function, /(x) 
1, in the interval (0, 1). 
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A. Mean Field Equation 



Let us consider a system of particles enclosed in a spherical region of radius R and volume 
V = A/3ttR^, with a total mass M distributed according to a smooth density p{r), normal- 
ized such that / d^rp{r) = 1, and interacting via a two body central potential (j){\r — r'\). 
If the potential is smooth, the entropy per particle in the microcanonical ensemble can be 
written in terms of the particle density as 

S = - J d^rp{r) [\nVp{r) - 1 ] + ^ In ( E - $ ) , (8) 

where E is the total energy and $ is the potential energy: 

= ^ y" d^rd^r' p{r)(t){\r -r'\) p{r') . (9) 

The volume V entering the first term of the r.h.s in has been included to make the entropy 
look dimensionally correct, and plays no significant role since it only shifts the entropy by a 
constant. The physical density is the absolute maximum of (§) under the constraint J p = 1. 
Differentiating with respect to p we arrive at the following integral equation: 

\nVp{r) = rfV0(|r-r'|)p(r'), (10) 

where f3 = 1/{E — ^) and p is the Lagrange multiplier for the constraint J p = 1. Defining 
u{r) by 

p{r) = ^ exp[/i + z/(r)] , (11) 



the constraint is solved by taking 



^'^ Jd^re'^C^)' ^^^^ 
Substituting (0) and (|12D in ([T0|), we obtain for z/(r): 

3^ /ciV0(|r-r-|)e-(^') 

= -^P Trfv^Kn • (13) 

If we take for 0(r) the Newtonian potential, we know that the entropy is not well defined. 
Nevertheless, it is still possible to start formally with the entropy functional which gives 
a finite result for any smooth distribution p(r), but is unbounded (see |1VC|) . There can 
still exist local maxima, which are then solutions of (|r^). By expanding the right-hand side 
of this equation in a series of spherical Bessel functions and truncating after N terms, one 
would obtain results equivalent to the ones we get using the softened potential. 

If we now particularize (|13D to the softened potential (^), we see that z/(r) obeys the 
same differential equation (^ as the potential. Imposing rotational symmetry on u we 
obtain the following ordinary differential equation: 
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= 0. (14) 

The general solution of this equation is a linear combination of {sm{ujkr)/r} and 
{cos(c(jfcr)/r}. The cosines should be absent from the solution since (pTSD implies that u 
is smooth. (Only at T = 0, i.e., (3 = oo, is u singular.) Indeed, it is shown explicitly in 
appendix A that the solution of (p!3D can be written as 



N 

^(r) = ^kMr/R), (15) 
fc=i 

where i^k are A'" numerical coefficients determined by the following set of equations: 

lyf^ = '^^'^^ -fo dxx^(f)k{x) exp{Efc^^fc0fc(x)} ^^g^ 
R Jo dxx'^ exp{Y.kJ^k(l)k{x)} 



The integral equation ([131) has been reduced to a system of non-linear algebraic equations 
with A^ unknowns. It can be solved by iteration, for instance with a Newton algorithm (see 
appendix B for a summary of the method used in this work). 



B. Thermodynamical quantities 



Using formula (^) of appendix A, it is straightforward to verify that, for a spherically 
symmetric mass distribution exp(/i + I'ir)), the potential energy is given by: 



R 



k=l 



/o" dr r2 



For an equilibrium distribution of the form ([TSD , Eq. ([Tq) implies 



1 R 



9/3^ GM^ 



N 

E 

k=l 



Hence, the total energy is 



E = - ^ ^ V 

P 9/?2 GM^ Y ' 

From (^ we easily obtain the equilibrium entropy: 
S = - \n(^J^ dxx"^ exp(^ Vk(t)k{x))^ + 



R 



GM^ 3/3 



-ln/3. 

2 ^ 



(17) 



(18) 



(19) 



(20) 



Since the entropy is stationary under variations of the mass distribution, the inverse tem- 
perature isl/T = dS/dE = [3 = 1/{E -^). 
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IV. RESULTS 



In order to present specific numerical results, it is convenient to work with dimensionless 
quantities. We measure the energy in units of the characteristic energy, GM'^/R, where 
M is the total mass and R the radius of the confining sphere. The dimensionless energy 
is then e = ER/ {GM'^). Any other quantity with dimensions of energy (like the tempera- 
ture 1/(3 and the potential energy $) must be also understood to be expressed in units of 
GM'^/R, and, similarly, magnitudes with dimensions of length are given in units of R. As a 
matter of terminology, we shall call Newtonian potential (NP) to the unregularized poten- 
tial, —GM'^/r, Newtonian entropy (NE) to its corresponding entropy, regularized potential 
(RP) to the potential regularized by Eq. (^), and regularized entropy (RE) to its associated 
entropy. 



A. N = 10 

For values of not too large, computations are easy. Let us describe the case = 10 
in detail. (See appendix B for a summary of the numerical methods used in this work). The 
solution of Eq. (0) as a function of the energy e provides all thermodynamic functions. For 
each value of e we found only one solution, which should then be the absolute maximum 
of (H). We shall return to this point later on, in Sec. [IV Q . Fig. ^ displays the inverse 
temperature (3 = 1/T versus e. In the thermodynamics of stable systems, this function must 



be monotonically decreasing, since the entropy is a convex function of the energy ||2^. In 
the present case, however, j3 decreases with e in the low and high energy regimes, but it 
increases for e in (—4.46, —0.2) and, consequently, the specific heat is negative in this energy 
interval. This is a consequence of the instability of the system. 

As is usually the case with these systems , the negative specific heat region is associated 
with a transition to a collapsed phase. To investigate this, let us define an order parameter 
K, = Ro/R, where Rq is the radius of the sphere centered at the origin which contains 95% of 
the mass (of course the value of 95% is arbitrary). Fig. ^ displays k versus e. At e = cxd the 
mass is distributed homogeneously, and then k = (0.95) ~ 0.9830. When the energy is 
reduced, k decreases monotonically and slowly. Notice the anomaly at e ~ —0.335; we shall 



discuss it in Sec. |IV B[ The collapsing order parameter n varies abruptly in the region where 
the specific heat is negative. It decays from k, ~ 0.95, corresponding to an homogeneous 
phase, to K ~ 0.1. In the later case, the mass distribution consists of a small dense core and 
an homogeneous tenuous halo. 

The results of this section are similar to those found by regularizing the potential with 
hard core repulsions and to those derived from the Lynden-Bell statistics applied to the 
unregularized potential |^|. It is remarkable that different regularizations lead to similar 
results. 



B. Newtonian potential 

It is interesting to compare the maximum of the RE with the local maximum of the NE. 
Substituting 0(|r — r'|) = — GM^/|r — r'\ in (|T3|), and using the fact that this potential is 
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a Green function of the Laplacian, we get the following differential equation 



V^u{r) = -^/5e^ exp(z/(r)), (21) 
which, for spherically symmetric u, is equivalent to the isothermal Lane-Emden equation 



^^ + H;?^ + |^,.e"« = (22) 
ar^ r dr 2 

The proper solutions of Eq. (0), with [3 and /i such that (3 = l/(e — $) and /i = 
— In/ drr^ exp z/(r), give local maxima of the entropy if e > —0.335 0,0]. The high energy 
phase should depend only weakly on the form of the potential at short distances. Therefore, 
the maximum of the RE in the high energy phase might be an approximation to the local 
maximum of the NE given by Eq. (p2|). This is indeed the case. 

To see how close the maximum of the RE is to the appropriate solution of eq. p^), we 
define a distance between functions hy D = max{|z/jv=io(?') — J^LE{r)\}, where the subscripts 
indicate the solutions of Eq. (|I6]) with = 10, and of Eq. (|22| ) respectively. For e > —0.335, 
i.e., when the Lane-Emden Eq. determines a local maximum of the NE, D < 10~^. The 
absolute maximum of the RE is indeed a very good approximation to the local maximum of 
the NE. 

Now, we can understand the anomaly in k around e ~ —0.335, which was mentioned 



in Sec. [IV A| and which can be appreciated in Fig. |^. At this point, which is close to the 
energy at which the solutions of the Lane-Emden eq. cease to be local maxima of the NE, 
the nature of the maximum of the RE also changes, originating anomalies such as the pick 
in 1/T (Fig. |) and the fissure in k (Fig. |]). 

The effect of the regularization is to deform the entropy functional dramatically for mass 
distributions p(r) which are very concentrated at the origin. These distributions get a huge 
amount of negative entropy after softening the potential, at least for A^ = 10, in such a way 
that there are no maxima of the RE close to them. On the other hand, the entropy of smooth 
distributions which are not concentrated is sensitive to the global form of the potential rather 
than to the short distance details. Therefore, these distributions have similar NE and RE, 
and they essentially do not feel the regularization. The solutions of the Lane-Emden Eq. 
belong to this class, and, consequently, close to them there is a local maximum of the RE 
which is indeed the global maximum for not too large A^, in particular for A^ = 10. 



C. A^ dependence 

The NP can be arbitrarily well approximated at short distances by a RP with A^ suf- 
ficiently large. Consequently, the maximum of the RE close to the local maximum of the 
NE will attain the later in the N oo limit, and, obviously, must cease to be the global 
maximum of the RE and become a local one for some value of A^, which will be denoted 
by A'c. Since the maximum of the entropy depends on the energy, A'c is a function of e. 
In principle, we can compute Nc{e) by solving Eq. (|l^) for large values of A^. In practice, 
however, this is very difficult and we must content ourselves with an estimate of Nc. 
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To get the estimate, let us first analyze how matter distributions with arbitrarily high 
NE can be built. There is an upper bound for the entropy functional (j^) if the potential 
energy (||) is bounded from below ($ > for any p(r)): 

S < 1 + ^ ln(e - . (23) 

In the case of the RP, $™^i„ = — A^. Since the entropy has an upper bound, it is reasonable 
to assume that it has a global maximum given by a regular function i/(r) of the form (|l5|), 
with coefficients z/^ verifying (p^G]). The potential energy associated to the NP has no lower 
bound and therefore = — oo. Hence, (p3D does not provide an upper bound for the NE. 
Indeed, it is straightforward to verify that the distribution 

f -r^ < r < ro 

^"■^ ^ { r„ < r<^ 

with < a < 1, has arbitrarily large entropy when we take the limit tq 0, while 
maintaining a In rg constant^. This is true for any value of e. We shall call these distributions, 
for any values of a and tq, special distributions (SD). If is large enough, there are SD 
with larger RE than the maximum close to the solution of the Lane-Emden equation.. 

As already claimed, it is very difficult to get the solutions of eq. ([T6|) for large values of 
A^. To overcome this problem and obtain an estimate of Nc, we shall study the restriction 
of the RE to particle distributions of the form (^) (SD). In such a way, we have a RE 
which depends only on two parameters, a and tq. Now, the maximization of this entropy 
with respect to a and vq is an easy task, even for very large values of A^. Obviously, the 
smaller value of A^ for which the maximum of the restricted RE is larger than the RE of 
the corresponding solution of the Lane-Emden equation will give an estimate of Nc- Strictly 
speaking, this estimate is an upper bound on Nc- 

Before analyzing the behavior with A^, let us look again to the A^ = 10 case, for which 
computations are easy. Fig. ^ displays the maximum of the restricted RE and the RE of the 
corresponding solution of ([T6|) , both for A^ = 10, as a function of e. The later distribution 
has always larger RE than any SD. This, besides the fact that we did not find other solutions 
by varying the initial guess, confirms that for each e only a local maximum of the A^ = 10 
RE exists. It is, obviously, the global maximum. 

To investigate the behavior with A^, we computed the estimate of the critical Nc for 
several values of e. As it could have been anticipated, grows with e. Table ^ shows the 
results. Column one displays e, column two the entropy of the solutions of the Lane-Emden 
eq., column three the maximum of the restriction of the RE to SD for the estimated Nc, and 



^Notice that in the limit rg ^ with alnro constant only an infinitesimal amount of matter 
collapses, while the rest is homogeneously distributed. This reflects the fact that it is enough 
that two particles (hard binaries) become arbitrarily close to make the potential energy arbitrarily 
negative, and therefore the kinetic energy arbitrarily large. However, as they constitute only two 
degrees of freedom, their contribution to the purely configurational entropy term / p(ln Vp — 1) of 
(^ is negligible. 
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column four the estimated Nc. It is apparent that, in the high energy phase (e > —0.335), 
we must go to N larger than 30 to see global maxima different from the solutions of the 
Lane-Emden equation. It is worth noting that a similar scenario was rigorously established 
by Kiessling for the equilibrium state of self-gravitating systems in contact with a thermal 
bath |]. 



V. CONCLUSIONS 

To define the thermodynamics of gravitational systems properly, the Newtonian potential 
must be regularized at short distances, removing its singularity. Only then is the entropy 
well defined, or, in a mean field approach, the entropy functional upper bounded. One way 
to introduce a regularization is by softening, i.e., by making the potential smooth at short 
distances while keeping it basically unchanged at long distances. There are infinitely many 
ways to achieve that. One interesting possibility is given by the truncation of the expansion 
of the gravitational potential in spherical Bessel functions to a given order A^, as in Eq. @. 
This regularization has the virtue of reducing the mean field integral equation to a system 
of algebraic equations with unknowns. This simplifies considerably the solution of the 
problem. 

The result which emerges from this approach is the following: if the regularization is 
mild enough, < 30, the system undergoes a phase transition separating a high energy 
homogeneous phase from a low energy collapsed phase. In the high energy phase, the mass 
distribution and the thermodynamic quantities are those of an isothermal sphere. Quan- 
titatively, they are very close to the solutions of the corresponding Lane-Emden equation. 
The low energy phase is characterized by a mass distribution consisting of a dense core 
surrounded by a tenuous halo. As usual in these cases 0, the transition from the HP to 
the CP takes place in an energy interval with negative specific heat, an indication of the 
thermodynamical instability of the system. These results are remarkably similar to those 
found with a different regularization (hard core spheres) and with those derived from 
the Lynden-Bell statistics applied to the unregularized potential p5| . We can then conclude 



that the thermodynamics is not very sensitive to the form of the regularization. 

The effect of a mild regularization is to deform the entropy functional in such a way that, 
in the high energy phase, the global maximum of the entropy with the regularized potential 
is very close to the local maximum with the unregularized potential. The analysis based on 
the Lane-Emden equation is therefore very accurate and the conclusions extracted from it 
hold. If the potential is too sharp at short distances, A^ > 30, we expect also a collapsing 
transition, which however will take place at a much higher energy than the one predicted 
by the analysis of the stability of the solutions of the Lane-Emden equation. Below the 
critical energy the global maximum of the RE will not be in the vicinity of the solution of 
the Lane-Emden equation, where, nevertheless, there will be a local maximum of the RE. 
Besides describing such metastable states, the Lane-Emden equation might be physically 
relevant in diluted self-gravitating systems with an interparticle distance high enough to be 
insensitive to a truncation of the expansion of the gravitational potential in spherical Bessel 
functions, Eq. (^, to 30 terms. 

Finally, let us comment on the structure of the low energy microcanonical equilibrium 
state when the regularized potential is very sharp (A^ > 30) at short distances. In this regime 



10 



there are high entropy mass distributions consisting of a small amount (infinitesimal when 
oo) of matter condensed and the rest homogeneously distributed. This might indicate 
that at these scales the system is not well described by a smooth density, and granularity is 
playing a major role. 



APPENDIX A: 



Let us show that the solutions of Eq. (|I3D for u rotationally symmetric are of the form 
(p!5|) if the potential is 0(r) = Y.k=i'Pk{r). Introducing spherical coordinates, eq. (|13D can 
be written 

3(3GM^ ^ dr' r' e"^"-'^ Jq dO sin ^^^(Vr^ + r'2 - 2rr' cos 6) 

""^'^ ~ h 2i>vvM • ^^^^ 

The integral in 6 can be readily performed, and gives 

,^ . ^ sin [cufc (r^ + r'^ — 2rr' cos 0)^/^1 „ , / n , / /n / , „n 

tUfc (r"^ + — 2rr' cosy) 

Eqs. (^) and imply 

'^('") — X]fc=i ^k'Pkix)-, with the coefficients z/^ determined by 

(0), QED. 



APPENDIX B: 

Since the numerical solution of Eq. (|1^) is central in this work, we shall outline in this 
appendix the method used to solve it. The problem is to find the roots of a vector function 
defined in a multidimensional space. Eq. ([16|) can be written as 

Fi{vi,...,VM) = (Bl) 

with i = 1, . . . ,N. We shall use matrix notation and denote by u the complete set 
(z/i, . . . , z/jv) and by F the vector (Fi, . . . , F^). 

To solve systems of equations like (|B1D we choose the Newton-Raphson method, which 



works as follows p9[: provided we have an initial guess u, which is close to the solution of 



we can expand Fj in Taylor series in a neighbourhood of v. 



N 

Fi{iy + 6iy) = Fi{u) + ^ Jij5uj + 0{5u^) (B2) 

i=i 

where Ja = is the Jacobian matrix. In matrix notation we have: 

F{u + 5v) = F{iy) + J-6u + 0{6u^) (B3) 

By neglecting terms of order higher than linear in 6u and by setting F{u + 6u) = 0, we 
obtain a set of linear equations for the corrections 6u that move each function Fi closer to 
zero simultaneously: 
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J.5v = -F (B4) 

This linear equation is a standar problem in numerical linear algebra and can be solved by 
LU decomposition. The corrections are then added to the initial guess, 

i^ne» = U + 6U (B5) 

and the process is iterated to convergence. It is possible to show that the method converges 
always provided the initial guess is close enough to the root. It can also spectacularly fail 
to converge indicating (though not proving) that the putative root does not exist nearby. 
To avoid problems with the poor global convergence of the method, we started with many 
different initial guess. We always found convergence to the same solution, except when 
was larger than 30, where we found only convergence at high energy. We never got two 
different solutions (whithin our convergence criterion, see below) by starting at two different 
points. 

Numerically, a convergence criterion is necessary. We stoped computations when one of 
this two conditions 

TV N 

J2\Sui\ < 10-^° or < 10"'° (B6) 

1=1 1=1 

was verified. 

Each time the function Fi was called, several integrals entering eq. ( [TB| ) were performed 



numerically, using a Romberg algorithm |29]. The integrands are smooth functions and it 



was possible to achieve a high precision with a relatively modest numerical effort. 
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TABLES 



Energy Entropy LE Entropy SD Nc 

Om -1.78525 -1.78247 79 

-0.12 -2.06341 -2.06124 56 

-0.20 -2.26269 -2.25417 44 

-0.30 -2.50812 -2.50519 32 
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FIG. 1. Newtonian and regularized potentials with iV = 10 and 20, in units of GM"^ jR. 
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FIG. 2. Inverse temperature 1/T versus energy e for AT = 10. 
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FIG. 3. Order parameter of collapsing (see text, Sec. [VA) versus energy e for = 10 
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FIG. 4. Entropy S versus energy e for A'^ = 10. The solid line is the entropy of the solution 
of Eq. ([l6|), and the dashed one correspond to the maximum entropy of distributions of the form 
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